Journal of Theoretical Biology

Trees accumulate somatic mutations throughout their long lifespan, resulting in genetic mosaicism among branches. While recent genomic studies quantified these mutations, they were largely limited to describing static patterns of variation. In this study, we developed a mathematical model to infer the dynamic processes of somatic mutation accumulation from snapshot genomic data obtained from four tropical trees (Dipterocarpaceae), which dominate tropical rain forests in Southeast Asia. Our model focus on genetic differences between shoot apical meristems (SAMs) at branch tips and explicitly incorporate stem cell dynamics within SAMs during shoot elongation and branching, enabling us to quantify somatic genetic drift arising from stem cell lineage replacement. By comparing model predictions with empirical data from Dipterocarpaceae trees, we estimated key parameters governing stem cell dynamics and somatic mutation rates. Our results indicate that both shoot elongation and branching involve replacement of stem cell lineages, leading to a moderate degree of somatic genetic drift. Accounting for stem cell dynamics resulted in slightly lower mutation rate estimates than previous approaches that ignored these processes. Using the estimated parameters, we further performed stochastic simulations to predict patterns of somatic mutations, including features not directly observed in the sampled trees, such as occasional deviations of somatic mutation phylogenies from physical architecture. Together, our modeling framework provides insights into how genetic mosaicism is shaped within tropical trees and reveals the stem cell dynamics underlying their long-term growth and accumulation of somatic mutations. (236 words) Highlights- We built mathematical models to predict the genetic differences between branch tips by somatic mutations. - The model considers the varying dynamics of stem cells in shoot meristem during shoot elongation and branching. - We compared the model prediction with empirical data from tropical trees, Dipterocarpaceae, and estimated the dynamics of stem cells and mutation rate. - Somatic mutation dynamics were shaped by somatic genetic drift arising from stem cell lineage replacement during shoot elongation and branching. - Accounting for stem cell dynamics led to slightly smaller estimates of mutation rates compared with previous estimates that ignored the dynamics. - Our models offer insights into how genetic variability is shaped in the tropical trees and the stem cell dynamics underlying their long-term growth.

Alzheimers disease (AD) is characterized by the accumulation of amyloid-{beta} (A{beta}), yet the specific link between plaque burden and cognitive decline remains a subject of intense investigation. This paper presents a mathematical model that simulates the coupled dynamics of A{beta} monomers, soluble oligomers, and fibrillar species in the brain tissue. By modifying existing moment equations to include a dedicated conservation equation for A{beta} monomers, the model explores how various microscopic processes, such as primary nucleation, surface-catalyzed secondary nucleation, fibril elongation, and fragmentation, contribute to macroscopic disease progression. Central to this study is the concept of "accumulated neurotoxicity" as a surrogate marker of biological age, defined as the time-integrated concentration of soluble A{beta} oligomers. Unlike plaque burden, accumulated neurotoxicity cannot be reversed, and the harm it causes depends critically on the sequence of events that produced it. Numerical results demonstrate that while plaque burden and neurotoxicity both increase over time, their relationship is non-linear and highly sensitive to the efficiency of protein degradation machinery. Specifically, impaired degradation leads to a rapid advancement of biological age relative to calendar age. The model further identifies oligomer dissociation and fibril fragmentation as potential protective mechanisms that can counterintuitively reduce neurotoxic burden by diverting monomers away from the soluble oligomer pool. These findings provide a quantitative framework for understanding why individuals with similar plaque burdens may experience vastly different cognitive outcomes, underscoring the importance of targeting soluble oligomers early in therapeutic interventions.

Patterns in the vegetation across arid and semiarid regions may be explained as a form of self-organization driven by water scarcity, and are often modeled through reaction-diffusion dynamics. Recent work has shown that similar mathematical models generate patterns on networks. However, these studies have focused on idealized topologies with no reference to natural pattern-forming systems. Our study aims at bridging these two fields: we employ a physical reaction-diffusion vegetation model, and gradually modify the topology of the diffusion network by adding random shortcuts over a 2-dimensional grid, interpolating between a regular lattice and a random network. We found that network topology strongly shapes both the resulting vegetation patterns and the precipitation range that supports them. Three behavioral regimes emerge. On a regular lattice, high-regularity patterns develop reflecting local diffusion processes. On a random network, the system is dominated by global pressure towards homogenization yielding either a uniform state or a single patch. In the intermediate shortcut density range, as the network topology resembles a small world network, the interaction between the two scales of diffusion generates two kinds of disordered patterns: low-regularity patterns with a well-defined characteristic wavelength, and irregular patterns characterized by a broad patch size distribution. These disordered patterns resemble real-world observations and, in our model, they show different responses to changing precipitation. Although we focused on dryland vegetation, we suggest that network-mediated diffusion could lead to similar mechanisms in a wide variety of pattern-forming systems. HighlightsO_LIWe study vegetation pattern formation over different diffusion network topologies. C_LIO_LITwo kinds of stable disordered patterns states develop over small world topologies. C_LIO_LILow-regularity patterns with a well-defined characteristic wavelength. C_LIO_LIIrregular patterns characterized by a broad patch size distribution. C_LIO_LIThese different kinds of disordered states show different relations to precipitation. C_LI

Mass-conserved reaction-diffusion systems are used as mathematical models for various phenomena such as cell polarity. Numerical simulations of this system present transient dynamics in which multiple stripe patterns converge to spatially monotonic patterns. Previous studies indicated that the transient dynamics are driven by a mass conservation law and by variations in the amount of substance contained in each pattern, which we refer to as "pattern flux". However, it is challenging to mathematically investigate these pattern dynamics. In this study, we introduce a reaction-diffusion compartment model to investigate the pattern dynamics in view of the conservation law and the pattern flux. This model is defined on multiple intervals (compartments), and diffusive couplings are imposed on each boundary of the compartments. Corresponding to the transient dynamics in the original system, we consider the dynamics around stripe patterns in the compartment model. We derive ordinary differential equations describing the pattern dynamics of the compartment model and analyze the existence and stability of equilibria for the reduced ODE with respect to the boundary parameters. For a specific parameter setting, we obtained results consistent with previous studies. Moreover, we present that the stripe patterns in the compartment model are potentially stabilized by changing the parameter, which is not observed in the original system. We expect that the methodology developed in this paper is extendable to various directions, such as membrane-induced pattern control.

The practical utility of many modern phylogenetic comparative methods can depend on how accurately mathematical models capture the evolutionary process of traits. Boucher and Demery (2016) described a new quantitative trait model, Brownian motion with reflective limits, that they anticipated might be of use in testing hypotheses about a particular sort of constraint on phenotypic character evolution. Since their analytic solution for the probability function under this bounded evolutionary scenario was not practical to evaluate for reasonably-sized trees, Boucher and Demery (2016) also identified a creative technique for computing the likelihood of their model. The basis of this methodology derives from the convergence of an equal-rates, symmetric, ordered Markov chain and continuous stochastic diffusion in the limit as the number of steps in our chain goes to {infty} (or, alternatively, as their widths decrease towards zero). We refer to this convergence in the limit as the discretized diffusion approximation or (more compactly) the discrete approximation. We realized that this discrete approximation of Boucher and Demery (2016) unlocked a number of additional models for the phylogenetic comparative analysis of discrete and continuous trait data, and we explore several of these in the present article. Specifically, we examine application of this discretized diffusion approximation to the threshold model from evolutionary quantitative genetics, to a new "semi-threshold" trait evolution model, to a joint model of discrete and continuous traits in which the discrete trait influences the rate of evolution of our continuous character, as well as a model where precisely the converse is true, and to a discrete character dependent multi-trend trended continuous trait evolution model. We conclude with some context for the origins of our article and discussion of other possible applications of this powerful approach.

Phenotypic diversity is often thought to arise from the evolutionary modification of developmental processes. However, developmental processes are tightly coupled in space and time, with each process beginning from conditions set by the one before it. While we know from dynamical systems theory that initial conditions can significantly affect a systems out-come, their importance as a source of phenotypic evolvability has been largely overlooked. Here we show for the first time, that phenotypic evolution can proceed through changes in developmental initial conditions while the underlying developmental process remains conserved. Somitogenesis is the process by which vertebral precursors, known as somites, are periodically patterned in the pre-somitic mesoderm (PSM). Somitic count (total number of somites) is thought to diversify through the evolution of components of somitogenesis such as the tempo of the segmentation clock or the mechanisms driving axial morphogenesis. Using two closely related species of Lake Malawi cichlid fishes that differ in vertebral counts, we show that somite count evolution has happened without changes to somitogenesis itself, but instead, by altering the size of the PSM at the onset of this process. This work will expand what we consider developmental drivers of phenotypic evolution and highlight the importance of comparative studies to understand the diversification of phenotypes.

We introduce a spectral existence criterion for the evolution of cooperation in the form of the inequality{lambda} maxb > c, where{lambda} max is the leading eigenvalue of an interaction operator encoding population structure, and b and c represent benefit and cost tradeoffs, respectively. Nowaks five rules for the evolution of cooperation correspond to cases in which the cooperation condition reduces to a scalar assortment coefficient. These results follow from the Price equation, which sheds light on a long-standing debate on the role of inclusive fitness and evolutionary dynamics in explaining the evolution of cooperation.

Human immunodeficiency virus (HIV) persistence remains a major barrier to cure due to the existence of long-lived latent reservoirs that evade immune clearance and persist despite combination antiretroviral therapy (ART). Although ART effectively suppresses viral replication, treatment interruption often leads to rapid viral rebound originating from these latent reservoirs. In this study, we develop a deterministic mathematical model describing the in vivo dynamics of HIV infection incorporating uninfected CD4+ T cells, infected cells, latent reservoirs, deep latent reservoirs, and infectious and non-infectious virions, while explicitly accounting for the therapeutic effects of reverse transcriptase inhibitors (RTIs), protease inhibitors (PIs), and Tat transcription inhibitors. Analytical results establish positivity and boundedness of solutions and derive the effective reproduction number Re using the next-generation matrix approach. Stability analysis shows that the virus-free equilibrium is locally asymptotically stable when Re < 1, while viral persistence occurs when Re > 1. Numerical simulations were performed to investigate therapy interactions, viral rebound following treatment interruption, and the impact of drug efficacy on viral set-points and latent reservoir dynamics. To further explore therapy interactions, three-dimensional viral set-point surfaces and heat maps were generated to examine how combinations of infection inhibition, viral production inhibition, and transcriptional inhibition influence viral dynamics. The simulations reveal that Tat inhibition suppresses viral transcription, thereby reducing the transition of infected cells into productive infection and limiting viral propagation when combined with conventional ART mechanisms. The therapy parameter planes further demonstrate that strong transcriptional inhibition promotes the transition of infected cells into deep latency, supporting the emerging block-and-lock strategy for functional HIV cure. In addition, a three-dimensional eradication boundary surface and therapy cube were constructed to identify regions of parameter space where Re < 1, corresponding to successful viral control. These visualizations show that viral eradication is unlikely when therapies act independently but becomes achievable when multiple therapeutic mechanisms act simultaneously. Overall, the results highlight the critical role of transcriptional inhibition through Tat-targeting therapies in complementing existing ART regimens. By simultaneously suppressing viral replication and promoting deep latency, Tat-based combination strategies may significantly reduce viral rebound and contribute to long-term functional control of HIV infection.

Most cancers are initiated by mutations that inactivate tumour-suppressor genes or activate oncogenes. Fitting a multistage model of carcinogenesis to the increase in cancer mortality with breed-specific size and lifespan in dogs has predicted that four somatic driver mutations are typically required to initiate cancer. This result is reconsidered here because it depends on the relationship between the number of at-risk cells and breed weight. Using a power function for this relationship results in higher quality models that support single driver mutations activating oncogenes. In addition, parameter estimates suggest that somatic mutation rates increase with weight, likely because of reduced investment in somatic maintenance. Regression of cancer mortality on body weight and lifespan shows that 56% of cancers in dogs are the result of mutations arising from somatic cell division, compared to 66% in humans. A further 7% of cancers may be due to inherited recessive mutations deactivating tumour suppressor genes, as indicated by the relationship between cancer mortality and breed germline homozygosity. Some of the remaining unexplained variation in cancer mortality may be explained by germline mutations underlying breed predispositions to specific cancers. Contrasting results with humans provide novel insights into the dynamics of carcinogenesis.

Cytokine-mediated communication is a central mechanism by which immune cells coordinate activation, differentiation and proliferation. While mechanistic reaction-diffusion models provide detailed descriptions of cytokine secretion and uptake at the cellular scale, their computational cost limits their applicability to large and densely packed cell populations. Previously employed approximations of cytokine diffusion fields rely on assumptions that neglect the influence of cellular geometry and volume exclusion. In this work, we study a macroscopic description of cytokine diffusion and reaction dynamics based on homogenization techniques, rigorously linking microscopic reaction-diffusion formulations to effective continuum models. The resulting homogenized equations replace discrete responder cells with a continuous density, while retaining essential features of cellular uptake and excluded-volume effects. Further, we show that in regimes with approximate radial symmetry, classical Yukawa-type solutions emerge as limiting cases of the homogenized model, provided appropriate correction factors are included. Overall, our approach allows efficient multiscale modeling of cytokine signaling in complex immune-cell environments.

Humans are just another species on Earth, but modern telecoupled societies and their socioeconomies impose immense consumption demands on the biosphere, detaching from common ecological rules. Starting from a simple ecological consumer-resource model, with humans as the consumers and terrestrial organic carbon (i.e., the biosphere) as the resource, we assume that technology modulates both human carrying capacity,{nu} 0, and the rate of biosphere consumption, 0. Three different functional-relation scenarios were tested, modulated by parameter a. In all three scenarios, equilibria and stability directly depended on the relative role that technology played in the model parameters, or the compound technological impact ({epsilon} {equiv} 0{nu}0). Moreover, two of the three scenarios showed Hopf bifurcations and regions with no equilibrium. The models were parameterized and fitted to actual data using a trajectory of more than 150 years. These analyses suggest that we are currently in a stable oscillatory spiral with no immediate Hopf bifurcation threat, but within a trajectory that continuously depletes the biosphere and approaches a collapse in human population size if no changes are made in the relationship that technology has with growth (i.e.,{nu} 0) versus consumption (i.e., 0) dynamics. Because our predatory dynamics also appear to have shifted from regular predator- prey dynamics toward a supply-demand scenario, with persistently increasing values, the threat of a Hopf bifurcation is now present in our trajectory: changes in the stability of the coexistence equilibrium may arise. This simple model warns that we must pay closer attention to the predatory relations that our technologies are creating with bio-sphere dynamics, in a way that goes beyond population numbers and technological development alone.

Qualitative models of Alzheimers pathology often posit that amyloid accumulation follows a sigmoid curve, indicating that the rate of deposition wanes over time. Longitudinal PET data now allow us to investigate amyloid accumulation trajectories with greater detail and over longer follow-up periods. We combine inferences from simulated amyloid trajectories, empirical PET data from the Alzheimers Disease Neuroimaging Initiative (ADNI), and the sampled iterative local approximation algorithm (SILA) to assess whether amyloid accumulation reaches a physiologic ceiling. We find that SILA reliably detects a ceiling, when present, across a range of simulated scenarios that impose a sigmoid shape. When fit to empirical data from ADNI, however, SILA does not appear to indicate the presence of a ceiling. Thus, we conclude that amyloid trajectories may not reach a physiologic ceiling during the stages of Alzheimers disease typically observed while patients remain under follow-up in cohort studies. Fits using SILA indicate that illustrative models of biomarker cascades, while useful tools for conceptualizing and interrogating pathologic processes, may not represent the shapes of amyloid trajectories accurately. Summary for General PublicAmyloid, a protein implicated in Alzheimers disease, is thought to reach a plateau in the brain, but methods that estimate how amyloid changes over time suggest it grows unabated. Gantenberg et al. use one such method and simulations to argue that amyloid does not reach a plateau during the typical course of Alzheimers.

Since most cancers are initiated by driver mutations arising in somatic cells, the risk of cancer should be explained by the probability of driver mutations as a function of the product of the number of stem cells and the stem cell division rate for a tissue. However, such models of driver-mutation initiated carcinogenesis fail to adequately predict cancer risk. It has been suggested that the missing component is greater adaptive cancer suppression in large tissues with high rates of stem cell division. This 8 may manifest as either a greater number of driver mutations required to initiate cancer or a lower stem cell mutation rate. These hypotheses are tested here using nonlinear regression to fit models that incorporate variation in the number of driver mutations and the stem cell mutation rate to data on the tissue-specific cancer risk, number of stem cells, and stem cell division rate. The greatest concordance between predicted and observed cancer risks is obtained by single driver mutations across tissues and stem cell mutation rates that decline with increasing lifetime numbers of stem cell divisions. This provides evidence of adaptive cancer suppression among tissues.

Interaction networks, in which nodes represent species and edges represent direct interactions between species, have a long and impactful history in community ecology. However, co-occurrence networks, where edges represent statistical relationships among species presences or abundances, are often easier to construct from lab and field data. It is clear that co-occurrence edges often do not represent direct interactions, but frameworks for the interpretation of co-occurrence networks have not kept pace with their generation. It is therefore unclear when and how these networks can be used to gain insight into community dynamics. Here, we use a Generalized Lotka-Volterra-based model to explore the contexts in which emergent properties of species interaction networks are identifiable in their resulting co-occurrence networks. We find that, in spite of many differences in direct edges, key features of the true interaction network, such as unipartite modularity, high-degree nodes (hubs), and bipartite modularity and nestedness, can be preserved in co-occurrence networks. In contrast, node degree distributions are not preserved even in the most idealized scenarios. We propose that networks derived from large co-occurrence datasets could therefore be used in future empirical work to test existing hypotheses of how emergent network structures drive ecological community dynamics.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.

Infectious disease surveillance systems in tropical countries show that respiratory disease incidence generally manifests as year-round activity with weak fluctuations and irregular seasonality. Previously, using a ten-year time series of influenza-like illness (ILI) collected from outpatient clinics in Ho Chi Minh City (HCMC), Vietnam, we found a combination of nonannual and annual signals driving these dynamics, but with unknown mechanisms. In this study, we use seven stochastic dynamical models incorporating humidity, temperature, and school term to investigate plausible mechanisms behind these annual and nonannual incidence trends. We use iterated filtering to fit the models and evaluate the models by comparing how well they replicate the combination of annual and nonannual signals. We find that a model including specific humidity, temperature, and school term best fits our observed data from HCMC and partially reproduces the irregular seasonality. The estimated effects from specific humidity and temperature on transmission are nonlinearly negative but weak. School dismissal is associated with decreased transmission, but also with low magnitude. Under these weak external drivers, we hypothesize that stochasticity makes a strong sub-annual cycle more likely to be observed in ILI disease dynamics. Our study shows a possible mechanism for respiratory disease dynamics in the tropics. When the external drivers are weak, the seasonality of respiratory disease dynamics is prone to the influence of stochasticity.

Evolutionary therapies regulate heterogeneous populations by altering selective pressures through treatment sequences in cancer and infections. This letter develops an invariant-set framework for treatment-induced containment based on positive triangular invariant sets. For periodically switched systems, sufficient conditions are derived for the existence of such invariant regions. Robustness with respect to mutation is established by showing that the invariant simplex persists under small perturbations of the subsystem matrices. In the two-phenotype case, the analysis yields an explicit mutation threshold that separates regimes in which therapy cycling maintains containment from regimes in which mutation can enable evolutionary escape. Simulations illustrate the geometry of the invariant sets and the role of mutation and dwell time in containment robustness.

We extend previous models of cumulative cultural evolution by incorporating structured populations and social networks. We examine how connectivity and network topology shape the accumulation of cultural complexity under unbiased (copy randomly), indirectly biased (copy successful individuals), and directly biased (copy successful traits) transmission. We consider random, scalefree, and small-world networks, as well as the communication structures introduced by Mason and Watts, and derive analytical approximations for the homogeneous case. We find that the effects of social structure depend strongly on the form of transmission bias. Under unbiased transmission, network effects are weak except at very low connectivity. Under indirect bias, cultural complexity increases with connectivity, whereas direct bias shows optimal performance at intermediate connectivity, reflecting a trade-off between diffusion and diversity. Differences across topologies are generally modest once the average degree is fixed. Overall, our results show that no single social structure universally promotes cumulative cultural evolution; instead, its effects depend primarily on the dynamics of learning and innovation.

Spatio-temporal interactions shape microbial community dynamics. Metabolism, through competition and cross-feeding, is a foundational mechanism of these interactions. Flux balance analysis enables efficient simulation of steady-state metabolism. Integrating these simulations through time, using dynamic flux balance analysis, provides temporal predictions of growth and metabolism. Incorporating spatial context, through partial differential equations, enables spatio-temporal simulation of microbial communities. In this chapter, we step through this sequential process, moving from steady-state, to temporal, to spatio-temporal simulation of microbial community metabolism. We provide an illustrative example using the modeling software COMETS (Computation of Microbial Ecosystems in Time and Space) to simulate interacting bacterial colonies of Bifidobacterium longum subsp. infantis and Anaerobutyricum hallii (previously Eubacterium hallii). Within this simulation, both competition and cross-feeding influenced the production of butyrate leading to an intermediate optimal interaction distance for metabolite production. We outline each step and provide open-source code such that this simulation can serve as a template for future spatio-temporal simulations of microbial community metabolism.

Tuberculosis (TB) remains a leading cause of death globally, with early mortality often driven by severe malnutrition and human immuno-deficiency virus (HIV) co-infection. Traditional survival analyses identify risk factors but remain associative, failing to capture the dynamic physiological collapse preceding death. In a novel interdisciplinary adaptation, we applied the Merton jump-diffusion structural framework from quantitative finance to model survival as a state of biological solvency, in which mortality occurs when a stochastic health trajectory crosses a critical failure threshold. We analysed a retrospective cohort of 15,182 TB patients in Cameroon over two decades. Adjusted body mass index (BMI) was conceptualized as a proxy for health capital and modeled using a stochastic process accounting for individual recovery trends, physiological instability, and acute clinical shocks. The study included predominantly young adult males (median age: 33 years) with a median BMI of 20.7 kg/m2. HIV co-infection was present in 35% of patients. The overall mortality rate during the 240 days follow-up period was 7.0%, with 55.1% of deaths occurring within the first 30 days. The model identified a critical failure threshold at BMI 17.329 kg/m2. HIV co-infection emerged as a key driver of metabolic instability, significantly increasing physiological volatility. Statistical validation confirmed that sudden clinical shocks were necessary to explain observed mortality patterns. The resulting Distance-to-Death (DtD) metric slightly outperformed standard associative extended Cox models in predicting survival, achieving a higher discriminative ability in testing set (Harrells C-index: 0.781 vs. 0.772; p = 0.038). Patients stratified into the highest-risk category showed a mortality rate of 16.7%, compared with 1.6% in the most stable group.This study bridges financial engineering and clinical epidemiology, offering a mechanistic understanding of how physiological reserves and metabolic instability determine survival. To support clinical application, we developed an interactive digital triage tool enabling identification of high-risk patients in resource-limited settings. Author summaryTuberculosis remains a major cause of death worldwide, particularly in people with poor nutrition or co-infection with HIV. In this study, we explored a new way to understand why some patients survive while others do not. We adapted a method originally used in finance to track the "health reserves" of patients over time, using body weight and related measures to estimate how close someone is to a critical health threshold. Our approach captures both gradual health decline and sudden medical complications, such as severe infections or rapid deterioration. By applying this method to a large group of patients in Cameroon, we found that a very low body weight is a strong warning sign for impending death and that HIV infection makes health outcomes less predictable. We also created a simple scoring tool that can help doctors identify patients at greatest risk, so that life-saving interventions and closer monitoring can be prioritized. This work bridges mathematical modeling and clinical care, offering a new way to assess patient vulnerability and improve outcomes in resource-limited settings.